model_fitting
load obs
obs_only <-
read_csv(glue("{params$model_dir_str}/data/stimulation_obvs.csv")) %>%
mutate(subj = as_factor(subj),
obs_degree = error,
error = obs_degree * (pi/180))## Parsed with column specification:
## cols(
## subj = col_double(),
## subj_index = col_double(),
## stimulation = col_double(),
## error = col_double()
## )# obs_only <- sim_data %>%
# unnest(subj_obs) %>%
# select(subj_index = subj,
# error = obs_radian,
# obs_degree,
# stimulation)peek at data
pre
obs_only %>%
filter(stimulation == 0) %>%
ggplot(aes(x = obs_degree)) +
geom_histogram(binwidth = 10, aes(y=..density..)) +
geom_rug() +
geom_density(aes(y=..density..)) +
facet_wrap(vars(subj), ncol = 1)
post
obs_only %>%
filter(stimulation == 1) %>%
ggplot(aes(x = obs_degree)) +
geom_histogram(binwidth = 10, aes(y=..density..)) +
geom_rug() +
geom_density(aes(y=..density..)) +
facet_wrap(vars(subj), ncol = 1)
fit brms
source(glue("{params$model_dir_str}/model_prior.R"))
print(bprior_full)## prior class coef group resp dpar nlpar bound
## 1 normal(3.8, 0.5) b intercept circSD
## 2 normal(0, 0.5) b stimulation circSD
## 3 normal(0, 0.5) sd Intercept subj circSD
## 4 normal(0, 0.5) sd stimulation subj circSD
## 5 normal(0, 1.5) b intercept theta
## 6 normal(0, 1.5) b stimulation theta
## 7 normal(0, 1) sd Intercept subj theta
## 8 normal(0, 1) sd stimulation subj thetaiter = 4000
warmup = 2000
cores = 4
chains = 4
n_post_samples = (iter - warmup) * chains
writeLines(
make_stancode(bform_full, obs_only, family = vm_uniform_mix, prior = bprior_full, stanvars = stanvars),
glue("{params$save_dir_str}/stan_code.txt")
)
model_fit <- brm(bform_full, obs_only, family = vm_uniform_mix, prior = bprior_full, stanvars = stanvars,
sample_prior = "yes",
warmup = warmup, iter = iter, cores = cores, chains = chains,
control = list(adapt_delta = 0.99), inits = 0,
file = glue("{params$save_dir_str}/obs_model_fit"))
print(model_fit)## Family: vm_uniform_mix
## Links: mu = identity; circSD = log; theta = logit; a = identity; b = identity
## Formula: error ~ 0
## circSD ~ 0 + intercept + stimulation + (1 + stimulation || subj)
## theta ~ 0 + intercept + stimulation + (1 + stimulation || subj)
## a = -3.14
## b = 3.14
## Data: obs_only (Number of observations: 504)
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup samples = 8000
##
## Group-Level Effects:
## ~subj (Number of levels: 2)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(circSD_Intercept) 0.36 0.25 0.02 0.96 1.00 3306
## sd(circSD_stimulation) 0.35 0.25 0.02 0.96 1.00 4334
## sd(theta_Intercept) 0.88 0.48 0.19 2.05 1.00 5352
## sd(theta_stimulation) 0.58 0.49 0.02 1.81 1.00 4306
## Tail_ESS
## sd(circSD_Intercept) 3035
## sd(circSD_stimulation) 3312
## sd(theta_Intercept) 4323
## sd(theta_stimulation) 3414
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## circSD_intercept 3.38 0.27 2.91 4.02 1.00 4535
## circSD_stimulation 0.26 0.29 -0.36 0.81 1.00 4995
## theta_intercept 0.10 0.62 -1.23 1.40 1.00 3782
## theta_stimulation -0.11 0.59 -1.31 1.14 1.00 4571
## Tail_ESS
## circSD_intercept 4303
## circSD_stimulation 4531
## theta_intercept 3996
## theta_stimulation 3529
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## b_circSD_intercept 3.38 0.27 2.91 4.02 1.00 4535
## b_circSD_stimulation 0.26 0.29 -0.36 0.81 1.00 4995
## b_theta_intercept 0.10 0.62 -1.23 1.40 1.00 3782
## b_theta_stimulation -0.11 0.59 -1.31 1.14 1.00 4571
## Tail_ESS
## b_circSD_intercept 4303
## b_circSD_stimulation 4531
## b_theta_intercept 3996
## b_theta_stimulation 3529
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).fit check
divergences
#check neff and rhat and divergences
np <- nuts_params(model_fit)
rhat <- brms::rhat(model_fit)
neff_rat <- neff_ratio(model_fit)
np %>%
filter(Parameter == "divergent__") %>%
summarise(n_div = sum(Value))## n_div
## 1 0rhat
mcmc_rhat(rhat) + yaxis_text(hjust = 1) + scale_x_continuous(breaks = pretty_breaks(6))## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
neff ratio
mcmc_neff(neff_rat) + yaxis_text(hjust = 1)
trace plots
mcmc_trace(as.array(model_fit$fit))
other
plot(model_fit, ask = FALSE)
* Notes
no divergences,
good ESS,
good rhat,
good fuzzy plots
plot posteriors
arrange samples
# compute summaries for plot
group_level_samples <-
spread_draws(model_fit, `(b|sd)_.*`, regex = TRUE) %>%
mutate(
# group level parameters
circSD_pre_mean = exp(b_circSD_intercept),
circSD_post_mean = exp(b_circSD_intercept + b_circSD_stimulation),
circSD_ES_mean = circSD_post_mean - circSD_pre_mean,
pMem_pre_mean = inv_logit(b_theta_intercept),
pMem_post_mean = inv_logit(b_theta_intercept + b_theta_stimulation),
pMem_ES_mean = pMem_post_mean - pMem_pre_mean,
# predicitve dist for group level parameters
circSD_pre_pred = exp(rnorm(n(), b_circSD_intercept, sd_subj__circSD_Intercept)),
circSD_post_pred = exp(rnorm(n(), b_circSD_intercept, sd_subj__circSD_Intercept) +
rnorm(n(), b_circSD_stimulation, sd_subj__circSD_stimulation)),
circSD_ES_pred = circSD_post_pred - circSD_pre_pred,
pMem_pre_pred = inv_logit(rnorm(n(), b_theta_intercept, sd_subj__theta_Intercept)),
pMem_post_pred = inv_logit(rnorm(n(), b_theta_intercept, sd_subj__theta_Intercept) +
rnorm(n(), b_theta_stimulation, sd_subj__theta_stimulation)),
pMem_ES_pred = pMem_post_pred - pMem_pre_pred
) %>%
select(-contains("b_"), -contains("sd_subj")) %>%
pivot_longer(-contains("."), names_to = c("param", "stat"), names_pattern = "(.*)_(.*)", values_to = "value") %>%
pivot_wider(names_from = stat, values_from = value)
group_level_summary <-
group_level_samples %>%
group_by(param) %>%
median_qi(.width = c(.5, .8, .95))
circSD_subj_samples <-
model_fit %>%
spread_draws(b_circSD_intercept, b_circSD_stimulation, r_subj__circSD[subj, term]) %>%
ungroup() %>%
pivot_wider(names_from = term, values_from = r_subj__circSD, names_prefix = "offset_") %>%
mutate(subj = subj,
circSD_pre = exp(b_circSD_intercept + offset_Intercept),
circSD_post = exp(b_circSD_intercept + offset_Intercept + b_circSD_stimulation + offset_stimulation),
circSD_ES = circSD_post - circSD_pre) %>%
select(-c(b_circSD_intercept, offset_Intercept, b_circSD_stimulation, offset_stimulation)) %>%
pivot_longer(contains("circSD"), names_to = "param", values_to = "value")
circSD_subj_summary <-
circSD_subj_samples %>%
group_by(subj, param) %>%
median_qi(.width = c(.90, .95))## Warning: unnest() has a new interface. See ?unnest for details.
## Try `df %>% unnest(c(.lower, .upper))`, with `mutate()` if needed
## Warning: unnest() has a new interface. See ?unnest for details.
## Try `df %>% unnest(c(.lower, .upper))`, with `mutate()` if neededpMem_subj_samples <-
model_fit %>%
spread_draws(b_theta_intercept, b_theta_stimulation, r_subj__theta[subj, term]) %>%
ungroup() %>%
pivot_wider(names_from = term, values_from = r_subj__theta, names_prefix = "offset_") %>%
mutate(subj = subj,
pMem_pre = inv_logit(b_theta_intercept + offset_Intercept),
pMem_post = inv_logit(b_theta_intercept + offset_Intercept + b_theta_stimulation + offset_stimulation),
pMem_ES = pMem_post - pMem_pre) %>%
select(-c(b_theta_intercept, offset_Intercept, b_theta_stimulation, offset_stimulation)) %>%
pivot_longer(contains("pMem"), names_to = "param", values_to = "value")
pMem_subj_summary <-
pMem_subj_samples %>%
group_by(subj, param) %>%
median_qi(.width = c(.90, .95))## Warning: unnest() has a new interface. See ?unnest for details.
## Try `df %>% unnest(c(.lower, .upper))`, with `mutate()` if needed
## Warning: unnest() has a new interface. See ?unnest for details.
## Try `df %>% unnest(c(.lower, .upper))`, with `mutate()` if neededgroup_level_samples %>%
select(-pred) %>%
group_by(param) %>%
median_qi(.width = c(.95))## Warning: unnest() has a new interface. See ?unnest for details.
## Try `df %>% unnest(c(.lower, .upper))`, with `mutate()` if needed## # A tibble: 6 x 7
## param mean .lower .upper .width .point .interval
## <chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
## 1 circSD_ES 8.88 -9.26 38.7 0.95 median qi
## 2 circSD_post 37.6 18.5 83.7 0.95 median qi
## 3 circSD_pre 28.7 18.3 55.7 0.95 median qi
## 4 pMem_ES -0.0282 -0.283 0.239 0.95 median qi
## 5 pMem_post 0.494 0.152 0.849 0.95 median qi
## 6 pMem_pre 0.523 0.226 0.802 0.95 median qigroup level posteriors
circSD_p1 <- group_level_samples %>%
filter(str_detect(param, "circSD")) %>%
ggplot() +
# posterior dist + interval for group mean
geom_halfeyeh(aes(y = param, x = mean), .width = c(.90, .95), position = position_nudge(y = 0.15)) +
# posterior predictive distribution for group means
stat_intervalh(aes(y = param, x = pred), .width = c(.5, .8, .95)) +
# posterior medians for each parameter estimate per subj
geom_point(data = circSD_subj_summary, aes(y = param, x = value), size = 2) +
# decorations
scale_color_brewer() +
scale_x_continuous(breaks = pretty_breaks(10)) +
coord_cartesian(xlim=c(-50, 150)) +
labs(subtitle = "circSD: group level mean posterior (median, 90%, 95% interval), \nsubject posterior medians, \ncondition predictive dist of subjects",
x = "circSD",
color = "interval")
# group level pMem pre, post and ES plot
pMem_p1 <- group_level_samples %>%
filter(str_detect(param, "pMem")) %>%
ggplot() +
# posterior dist + interval for group mean
geom_halfeyeh(aes(y = param, x = mean), .width = c(.90, .95), position = position_nudge(y = 0.15)) +
# posterior predictive distribution for group means
stat_intervalh(aes(y = param, x = pred), .width = c(.5, .8, .95)) +
# posterior medians for each parameter estimate per subj
geom_point(data = pMem_subj_summary, aes(y = param, x = value), size = 2) +
# decorations
scale_color_brewer() +
scale_x_continuous(breaks = pretty_breaks(10)) +
#coord_cartesian(xlim=c(-50, 150)) +
labs(subtitle = "pMem: group level mean posterior (median, 90%, 95% interval), \nsubject posterior medians, \ncondition predictive dist of subjects",
x = "pMem",
color = "interval")
plot_grid(circSD_p1, pMem_p1, align = "hv", ncol = 1)
circSD pre, post, ES
# circSD pre: group level posteriors and subject posteriors
circSD_p2 <-
ggplot() +
# plot group mean circSD_pre posterior and subject circSD_pre posterior
geom_halfeyeh(data = rbind(
group_level_samples %>%
filter(str_detect(param, "circSD_pre")) %>%
select(-pred, value = mean)
,
circSD_subj_samples %>%
filter(str_detect(param, "circSD_pre")) %>%
unite(param, param, subj) )
, aes(y = param, x = value), .width = c(.90, .95)) +
# plot pre condition group predictive distribution
stat_intervalh(data = group_level_samples %>% filter(str_detect(param, "circSD_pre")),
aes(y = param , x = pred),
.width = c(.5, .8, .95),
position = position_nudge(y = -0.15)
) +
# show modeled subject circSD_pre posterior medians in the prediction band
geom_point(data = circSD_subj_summary %>% filter(param == "circSD_pre"),
aes(y = param, x = value),
size = 2,
position = position_nudge(y = -0.15)) +
scale_color_brewer() +
scale_x_continuous(breaks = pretty_breaks(10)) +
labs(subtitle = "circSD_pre: group level mean posterior (median, 90%, 95% interval), \nsubject posterior, \ncondition predictive dist of subjects",
x = "circSD",
color = "interval")
circSD_p3 <-
ggplot() +
# plot group mean circSD_post posterior and subject circSD_post posterior
geom_halfeyeh(data = rbind(
group_level_samples %>%
filter(str_detect(param, "circSD_post")) %>%
select(-pred, value = mean)
,
circSD_subj_samples %>%
filter(str_detect(param, "circSD_post")) %>%
unite(param, param, subj) )
, aes(y = param, x = value), .width = c(.90, .95)) +
# plot post condition group predictive distribution
stat_intervalh(data = group_level_samples %>% filter(str_detect(param, "circSD_post")),
aes(y = param , x = pred),
.width = c(.5, .8, .95),
position = position_nudge(y = -0.15)
) +
# show modeled subject circSD_post posterior medians in the prediction band
geom_point(data = circSD_subj_summary %>% filter(param == "circSD_post"),
aes(y = param, x = value),
size = 2,
position = position_nudge(y = -0.15)) +
scale_color_brewer() +
scale_x_continuous(breaks = pretty_breaks(10)) +
labs(subtitle = "circSD_post: group level mean posterior (median, 90%, 95% interval), \nsubject posterior, \ncondition predictive dist of subjects",
x = "circSD",
color = "interval")
circSD_p4 <-
ggplot() +
# plot group mean circSD_ES posterior and subject circSD_ES posterior
geom_halfeyeh(data = rbind(
group_level_samples %>%
filter(str_detect(param, "circSD_ES")) %>%
select(-pred, value = mean)
,
circSD_subj_samples %>%
filter(str_detect(param, "circSD_ES")) %>%
unite(param, param, subj) )
, aes(y = param, x = value), .width = c(.90, .95)) +
# plot ES group predictive distribution
stat_intervalh(data = group_level_samples %>% filter(str_detect(param, "circSD_ES")),
aes(y = param , x = pred),
.width = c(.5, .8, .95),
position = position_nudge(y = -0.15)
) +
# show modeled subject circSD_ES posterior medians in the prediction band
geom_point(data = circSD_subj_summary %>% filter(param == "circSD_ES"),
aes(y = param, x = value),
size = 2,
position = position_nudge(y = -0.15)) +
scale_color_brewer() +
scale_x_continuous(breaks = pretty_breaks(10)) +
labs(subtitle = "circSD_ES: group level mean posterior (median, 90%, 95% interval), \nsubject posterior, \nES predictive dist of subjects",
x = "detla circSD",
color = "interval")
plot_grid(circSD_p2, circSD_p3, circSD_p4, ncol = 1, align = "hv")
pMem pre, post, ES
# pMem_ pre: group level posteriors and subject posteriors
pMem_p2 <-
ggplot() +
# plot group mean pMem__pre posterior and subject cpMem__pre posterior
geom_halfeyeh(data = rbind(
group_level_samples %>%
filter(str_detect(param, "pMem_pre")) %>%
select(-pred, value = mean)
,
pMem_subj_samples %>%
filter(str_detect(param, "pMem_pre")) %>%
unite(param, param, subj) )
, aes(y = param, x = value), .width = c(.90, .95)) +
# plot pre condition group predictive distribution
stat_intervalh(data = group_level_samples %>% filter(str_detect(param, "pMem_pre")),
aes(y = param , x = pred),
.width = c(.5, .8, .95),
position = position_nudge(y = -0.15)
) +
# show modeled subject pMem_pre posterior medians in the prediction band
geom_point(data = pMem_subj_summary %>% filter(param == "pMem_pre"),
aes(y = param, x = value),
size = 2,
position = position_nudge(y = -0.15)) +
scale_color_brewer() +
scale_x_continuous(breaks = pretty_breaks(10)) +
labs(subtitle = "pMem_pre: group level mean posterior (median, 90%, 95% interval), \nsubject posterior, \ncondition predictive dist of subjects",
x = "pMem",
color = "interval")
pMem_p3 <-
ggplot() +
# plot group mean pMem_post posterior and subject pMem_post posterior
geom_halfeyeh(data = rbind(
group_level_samples %>%
filter(str_detect(param, "pMem_post")) %>%
select(-pred, value = mean)
,
pMem_subj_samples %>%
filter(str_detect(param, "pMem_post")) %>%
unite(param, param, subj) )
, aes(y = param, x = value), .width = c(.90, .95)) +
# plot post condition group predictive distribution
stat_intervalh(data = group_level_samples %>% filter(str_detect(param, "pMem_post")),
aes(y = param , x = pred),
.width = c(.5, .8, .95),
position = position_nudge(y = -0.15)
) +
# show modeled subject pMem_post posterior medians in the prediction band
geom_point(data = pMem_subj_summary %>% filter(param == "pMem_post"),
aes(y = param, x = value),
size = 2,
position = position_nudge(y = -0.15)) +
scale_color_brewer() +
scale_x_continuous(breaks = pretty_breaks(10)) +
labs(subtitle = "pMem_post: group level mean posterior (median, 90%, 95% interval), \nsubject posterior, \ncondition predictive dist of subjects",
x = "pMem",
color = "interval")
pMem_p4 <-
ggplot() +
# plot group mean pMem_ES posterior and subject pMem_ES posterior
geom_halfeyeh(data = rbind(
group_level_samples %>%
filter(str_detect(param, "pMem_ES")) %>%
select(-pred, value = mean)
,
pMem_subj_samples %>%
filter(str_detect(param, "pMem_ES")) %>%
unite(param, param, subj) )
, aes(y = param, x = value), .width = c(.90, .95)) +
# plot ES group predictive distribution
stat_intervalh(data = group_level_samples %>% filter(str_detect(param, "pMem_ES")),
aes(y = param , x = pred),
.width = c(.5, .8, .95),
position = position_nudge(y = -0.15)
) +
# show modeled subject pMem_ES posterior medians in the prediction band
geom_point(data = pMem_subj_summary %>% filter(param == "pMem_ES"),
aes(y = param, x = value),
size = 2,
position = position_nudge(y = -0.15)) +
scale_color_brewer() +
scale_x_continuous(breaks = pretty_breaks(10)) +
labs(subtitle = "pMem_ES: group level mean posterior (median, 90%, 95% interval), \nsubject posterior, \nES predictive dist of subjects",
x = "detla pMem",
color = "interval")plot_grid(pMem_p2, pMem_p3, pMem_p4, ncol = 1, align = "hv")
group level joint posteriors
group_level_samples %>%
pivot_wider(id_cols = contains("."), names_from = param, values_from = mean) %>%
select(-contains(".")) %>%
mcmc_pairs(off_diag_fun = "hex")## Warning: Only one chain in 'x'. This plot is more useful with multiple
## chains.
plot posteriors w/ priors
# compute summaries for plot with priors
group_level_samples <-
spread_draws(model_fit, `(prior_)?(b|sd)_.*`, regex = TRUE) %>%
mutate(
# prior
# group level parameters
prior_circSD_pre_mean = exp(prior_b_circSD_intercept),
prior_circSD_post_mean = exp(prior_b_circSD_intercept + prior_b_circSD_stimulation),
prior_circSD_ES_mean = prior_circSD_post_mean - prior_circSD_pre_mean,
prior_pMem_pre_mean = inv_logit(prior_b_theta_intercept),
prior_pMem_post_mean = inv_logit(prior_b_theta_intercept + prior_b_theta_stimulation),
prior_pMem_ES_mean = prior_pMem_post_mean - prior_pMem_pre_mean,
# predicitve dist for group level parameters
prior_circSD_pre_pred = exp(rnorm(n(), prior_b_circSD_intercept, prior_sd_subj__circSD_Intercept)),
prior_circSD_post_pred = exp(rnorm(n(), prior_b_circSD_intercept, prior_sd_subj__circSD_Intercept) +
rnorm(n(), prior_b_circSD_stimulation, prior_sd_subj__circSD_stimulation)),
prior_circSD_ES_pred = prior_circSD_post_pred - prior_circSD_pre_pred,
prior_pMem_pre_pred = inv_logit(rnorm(n(), prior_b_theta_intercept, prior_sd_subj__theta_Intercept)),
prior_pMem_post_pred = inv_logit(rnorm(n(), prior_b_theta_intercept, prior_sd_subj__theta_Intercept) +
rnorm(n(), prior_b_theta_stimulation, prior_sd_subj__theta_stimulation)),
prior_pMem_ES_pred = prior_pMem_post_pred - prior_pMem_pre_pred,
# posteriors
# group level parameters
circSD_pre_mean = exp(b_circSD_intercept),
circSD_post_mean = exp(b_circSD_intercept + b_circSD_stimulation),
circSD_ES_mean = circSD_post_mean - circSD_pre_mean,
pMem_pre_mean = inv_logit(b_theta_intercept),
pMem_post_mean = inv_logit(b_theta_intercept + b_theta_stimulation),
pMem_ES_mean = pMem_post_mean - pMem_pre_mean,
# predicitve dist for group level parameters
circSD_pre_pred = exp(rnorm(n(), b_circSD_intercept, sd_subj__circSD_Intercept)),
circSD_post_pred = exp(rnorm(n(), b_circSD_intercept, sd_subj__circSD_Intercept) +
rnorm(n(), b_circSD_stimulation, sd_subj__circSD_stimulation)),
circSD_ES_pred = circSD_post_pred - circSD_pre_pred,
pMem_pre_pred = inv_logit(rnorm(n(), b_theta_intercept, sd_subj__theta_Intercept)),
pMem_post_pred = inv_logit(rnorm(n(), b_theta_intercept, sd_subj__theta_Intercept) +
rnorm(n(), b_theta_stimulation, sd_subj__theta_stimulation)),
pMem_ES_pred = pMem_post_pred - pMem_pre_pred
) %>%
select(-contains("b_"), -contains("sd_subj")) %>%
pivot_longer(-contains("."), names_to = c( "param", "stat"), names_pattern = "(.*)_(.*)", values_to = "value") %>%
pivot_wider(names_from = stat, values_from = value)group level posteriors
circSD_p1_wPrior <-
group_level_samples %>%
filter(str_detect(param, "circSD")) %>%
mutate(param = fct_relevel(param, c("prior_circSD_ES", "circSD_ES", "prior_circSD_pre", "circSD_pre", "prior_circSD_post", "circSD_post"))) %>%
ggplot() +
# posterior dist + interval for group mean
geom_halfeyeh(aes(y = param, x = mean), .width = c(.90, .95), position = position_nudge(y = 0.15)) +
# posterior predictive distribution for group means
stat_intervalh(aes(y = param, x = pred), .width = c(.5, .8, .95)) +
# posterior medians for each parameter estimate per subj
geom_point(data = circSD_subj_summary, aes(y = param, x = value), size = 2) +
# prior dist + interval for group mean
#geom_halfeyeh(data = group_level_prior_param_samples, aes(y = param, x = values), .width = c(.90, .95)) +
# decorations
scale_color_brewer() +
scale_x_continuous(breaks = pretty_breaks(10)) +
coord_cartesian(xlim=c(-50, 150)) +
labs(subtitle = "circSD: group level mean prior + posterior (median, 90%, 95% interval), \nsubject posterior medians, \ncondition predictive dist of subjects",
x = "circSD",
color = "interval")
# group level pMem pre, post and ES plot
pMem_p1_wPrior <-
group_level_samples %>%
filter(str_detect(param, "pMem")) %>%
mutate(param = fct_relevel(param, c("prior_pMem_ES", "pMem_ES", "prior_pMem_pre", "pMem_pre", "prior_pMem_post", "pMem_post"))) %>%
ggplot() +
# posterior dist + interval for group mean
geom_halfeyeh(aes(y = param, x = mean), .width = c(.90, .95), position = position_nudge(y = 0.15)) +
# posterior predictive distribution for group means
stat_intervalh(aes(y = param, x = pred), .width = c(.5, .8, .95)) +
# posterior medians for each parameter estimate per subj
geom_point(data = pMem_subj_summary, aes(y = param, x = value), size = 2) +
# decorations
scale_color_brewer() +
scale_x_continuous(breaks = pretty_breaks(10)) +
#coord_cartesian(xlim=c(-50, 150)) +
labs(subtitle = "pMem: group level mean prior + posterior (median, 90%, 95% interval), \nsubject posterior medians, \ncondition predictive dist of subjects",
x = "pMem",
color = "interval")
plot_grid(circSD_p1_wPrior, pMem_p1_wPrior, align = "hv", ncol = 1)
* Notes
Main conclusion is that both group mean and subject predicitive dist show reductions in uncertainty, though lots of probability still over zero for both pMem and circSD ES.
There is a sign of a circSD increase, more data will make that clearer.
posterior predictive plot of errors
sim
conditions <- c(0,1)
nobs_per_cond_sim <- 1500
post_pred_fpath <- glue("{params$save_dir_str}/post_pred_obs.rds")
if (file.exists(post_pred_fpath)){
post_pred_sim <- readRDS(post_pred_fpath)
}else{
print("simulating")
posterior_arranged <-
spread_draws(model_fit, `(b|sd)_.*`, regex = TRUE) %>%
sample_n(2e3) %>%
mutate(sim_num = 1:n(),
nsubj = 1,
nobs_per_cond = nobs_per_cond_sim) %>%
select(-contains(".")) %>%
select(sim_num,
alpha0_mu = b_circSD_intercept,
alpha0_sigma = sd_subj__circSD_Intercept,
alphaD_mu = b_circSD_stimulation,
alphaD_sigma = sd_subj__circSD_stimulation,
beta0_mu = b_theta_intercept,
beta0_sigma = sd_subj__theta_Intercept,
betaD_mu = b_theta_stimulation,
betaD_sigma = sd_subj__theta_stimulation,
nsubj,
nobs_per_cond)
post_pred_sim <-
posterior_arranged %>%
mutate(
# use draw_subj to sample nsubj_sim per sim using group-level parameter draws
dataset = pmap(posterior_arranged, draw_subj),
stimulation = list(stimulation = rep(conditions, each = nobs_per_cond_sim))) %>%
# first unnest dataset, expanding by nsubj_sim and copying stimulation list to each subj
unnest(dataset) %>%
# then unnest stimulation, expanding by nobs_per_cond_sim*2
unnest(stimulation) %>%
# now use likelihood to simulation observations
mutate(
# evaluate and delink linear model on pMem
pMem = inv_logit(subj_beta0 + (subj_betaD * stimulation)),
# evaluate and delink linear model on circSD/kappa
k = sd2k_vec(
pracma::deg2rad(
exp(subj_alpha0 + (subj_alphaD * stimulation)))),
# use pMem to draw a 1 or 0 for each trial
memFlip = rbernoulli(n(), pMem),
# use k to draw from vonMises for each trial
vm_draw = rvonmises_vec(1, pi, k) - pi,
# draw from unif for each trial
unif_draw = runif(n(), -pi, pi),
# assign either vm_draw or unif_draw to each trial, depending on memFlip
obs_radian = memFlip * vm_draw + (1 - memFlip) * unif_draw,
# convert to degrees
obs_degree = obs_radian * (180/pi)
) %>%
select(-c(pMem, k, memFlip, vm_draw, unif_draw)) %>%
nest(subj_obs = c(stimulation, obs_degree, obs_radian)) %>%
nest(dataset = c(subj, subj_alpha0, subj_alphaD, subj_beta0, subj_betaD, nobs_per_condition, subj_obs))
saveRDS(post_pred_sim, post_pred_fpath)
}conditions <- c(0,1)
nobs_per_cond_sim <- 1500
post_predmean_fpath <- glue("{params$save_dir_str}/post_predmean_obs.rds")
if (file.exists(post_predmean_fpath)){
post_predmean_sim <- readRDS(post_predmean_fpath)
}else{
posterior_arranged <-
spread_draws(model_fit, `(b|sd)_.*`, regex = TRUE) %>%
sample_n(2e3) %>%
mutate(sim_num = 1:n(),
nobs_per_cond = nobs_per_cond_sim ) %>%
select(-contains(".")) %>%
select(sim_num,
subj_alpha0 = b_circSD_intercept,
subj_alphaD = b_circSD_stimulation,
subj_beta0 = b_theta_intercept,
subj_betaD = b_theta_stimulation,
nobs_per_cond)
post_predmean_sim <-
posterior_arranged %>%
# add stimulation covariates
mutate(
stimulation = list(stimulation = rep(conditions, each = nobs_per_cond_sim))
) %>%
# then unnest stimulation, expanding by nobs_per_cond_sim*2
unnest(stimulation) %>%
# now use likelihood to simulation observations
mutate(
# evaluate and delink linear model on pMem
pMem = inv_logit(subj_beta0 + (subj_betaD * stimulation)),
# evaluate and delink linear model on circSD/kappa
k = sd2k_vec(
pracma::deg2rad(
exp(subj_alpha0 + (subj_alphaD * stimulation)))),
# use pMem to draw a 1 or 0 for each trial
memFlip = rbernoulli(n(), pMem),
# use k to draw from vonMises for each trial
vm_draw = rvonmises_vec(1, pi, k) - pi,
# draw from unif for each trial
unif_draw = runif(n(), -pi, pi),
# assign either vm_draw or unif_draw to each trial, depending on memFlip
obs_radian = memFlip * vm_draw + (1 - memFlip) * unif_draw,
# convert to degrees
obs_degree = obs_radian * (180/pi)
) %>%
select(-c(pMem, k, memFlip, vm_draw, unif_draw)) %>%
nest(dataset = c(stimulation, obs_degree, obs_radian))
saveRDS(post_predmean_sim, post_predmean_fpath)
}posterior data prediction, w/ group-level variance
#######################################################
# calculate histogram quantile mats from each condition
sim_subj_obs_hist_count <- function(dataset, condition = 0){
dataset_obs <- dataset %>%
sample_n(1) %>%
unnest(subj_obs) %>%
filter(stimulation == condition) %>%
select(obs_degree)
breaks <- seq(-180, 180, 5)
bincount <- hist(dataset_obs$obs_degree, breaks = breaks, plot = FALSE)$counts
bincount_names <- glue("c{breaks[-1]}")
names(bincount) <- bincount_names
bincount_df <- data.frame(as.list(bincount))
return(bincount_df)
}
make_quantmat <- function(sim_datasets, condition = 0){
bincounts <- sim_datasets %>%
select(dataset) %>%
mutate(subj_hist_counts = map(dataset, sim_subj_obs_hist_count, condition)) %>%
select(-dataset) %>%
unnest(subj_hist_counts) %>%
as_tibble()
xvals <- seq(-177.5, 177.5, 5)
probs <- seq(0.1,0.9,0.1)
quantmat <- as.data.frame(matrix(NA, nrow=ncol(bincounts), ncol=length(probs)))
names(quantmat) <- paste0("p",probs)
quantmat <- cbind(quantmat, xvals)
for (iQuant in 1:length(probs)){
quantmat[,paste0("p",probs[iQuant])] <- as.numeric(summarise_all(bincounts, ~quantile(., probs[iQuant])))
}
return(quantmat)
}
quantmat_cond0 <- make_quantmat(post_pred_sim, 0)
quantmat_cond1 <- make_quantmat(post_pred_sim, 1)
#######################################################
# calculate ecdf quantile mats from each condition
# unnest post_pred_sim, using only 1 subj/dataset
unnested <-
post_pred_sim %>%
unnest(dataset) %>%
group_by(sim_num) %>%
sample_n(1) %>%
ungroup() %>%
unnest(subj_obs)
# calc quantiles mat for pre condition
ecdf_res_stim0 <-
unnested %>%
filter(stimulation == 0) %>%
group_by(sim_num) %>%
group_map(~ecdf(.$obs_degree )(seq(-180, 180, 1)))
stim0_ecdf_quantiles <- bind_cols(
tibble(x_val = seq(-180, 180, 1)),
as_tibble(colQuantiles(do.call(rbind, ecdf_res_stim0), probs = c(0.95, 0.5, 0.05 )))
)
# calc quantiles mat for post condition
ecdf_res_stim1 <-
unnested %>%
filter(stimulation == 1) %>%
group_by(sim_num) %>%
group_map(~ecdf(.$obs_degree )(seq(-180, 180, 1)))
stim1_ecdf_quantiles <- bind_cols(
tibble(x_val = seq(-180, 180, 1)),
as_tibble(colQuantiles(do.call(rbind, ecdf_res_stim1), probs = c(0.95, 0.5, 0.05 )))
)# blues
b_light <- "#8C9BC4"
b_light_highlight <- "#A0ADCE"
b_mid <- "#546BA9"
b_mid_highlight <- "#7385B8"
b_dark <- "#002381"
b_dark_highlight <- "#2E4B97"
#betancourt reds
r_light <- "#DCBCBC"
r_light_highlight <- "#C79999"
r_mid <- "#B97C7C"
r_mid_highlight <- "#A25050"
r_dark <- "#8F2727"
r_dark_highlight <- "#7C0000"
#######################################################
# plot histogram(density) per condition
ggplot(quantmat_cond0, aes(x = xvals)) +
geom_ribbon(aes(ymax = p0.9, ymin = p0.1), fill = r_light, alpha = 0.4) +
geom_line(aes(y = p0.5), color = r_dark, size = 1) +
geom_ribbon(data = quantmat_cond1, aes(ymax = p0.9, ymin = p0.1), fill = b_light, alpha = 0.4) +
geom_line(data = quantmat_cond1, aes(y = p0.5), color = b_dark, size = 1) +
scale_x_continuous(breaks=pretty_breaks(10)) +
labs(x = "error (degrees) [red = pre, blue = post]",
y = "count +/- quantile",
subtitle = glue("per-condition posterior pred dist (median line, 90% interval over {nrow(post_pred_sim)} sim datasets)\n({nobs_per_cond_sim} samples/cond, per subj-level draw, per group-level mean + sd draw)")) +
theme_cowplot()
#######################################################
# plot ecdf per condition
ggplot() +
geom_ribbon(data = stim0_ecdf_quantiles, aes(x = x_val, ymax = `95%`, ymin = `5%`), fill = "red", alpha = 0.3) +
geom_line(data = stim0_ecdf_quantiles, aes(x = x_val, y = `50%`), color = "red", size = 1) +
geom_ribbon(data = stim1_ecdf_quantiles, aes(x = x_val, ymax = `95%`, ymin = `5%`), fill = "blue", alpha = 0.3) +
geom_line(data = stim1_ecdf_quantiles, aes(x = x_val, y = `50%`), color = "blue", size = 1) +
scale_x_continuous(breaks=pretty_breaks(10)) +
geom_hline(yintercept = seq(0, 1, 0.25), linetype = "dashed", alpha = 0.2) +
labs(x = "error (degrees) [red = pre, blue = post]",
y = "cumulative prob.",
subtitle = glue("per-condition posterior pred cdf (median line, 90% interval over {nrow(post_pred_sim)} sim datasets) \n({nobs_per_cond_sim} samples/cond, per subj-level draw, per group-level mean + sd draw")) +
theme_cowplot()
posterior data prediction, w/o group-level variance
#######################################################
# calculate histogram quantile mats from each condition
sim_single_subj_obs_hist_count <- function(dataset, condition = 0){
dataset_obs <- dataset %>%
filter(stimulation == condition) %>%
select(obs_degree)
breaks <- seq(-180, 180, 5)
bincount <- hist(dataset_obs$obs_degree, breaks = breaks, plot = FALSE)$counts
bincount_names <- glue("c{breaks[-1]}")
names(bincount) <- bincount_names
bincount_df <- data.frame(as.list(bincount))
return(bincount_df)
}
make_quantmat <- function(sim_datasets, condition = 0){
bincounts <- sim_datasets %>%
select(dataset) %>%
mutate(subj_hist_counts = map(dataset, sim_single_subj_obs_hist_count, condition)) %>%
select(-dataset) %>%
unnest(subj_hist_counts) %>%
as_tibble()
xvals <- seq(-177.5, 177.5, 5)
probs <- seq(0.1,0.9,0.1)
quantmat <- as.data.frame(matrix(NA, nrow=ncol(bincounts), ncol=length(probs)))
names(quantmat) <- paste0("p",probs)
quantmat <- cbind(quantmat, xvals)
for (iQuant in 1:length(probs)){
quantmat[,paste0("p",probs[iQuant])] <- as.numeric(summarise_all(bincounts, ~quantile(., probs[iQuant])))
}
return(quantmat)
}
quantmat_cond0 <- make_quantmat(post_predmean_sim, 0)
quantmat_cond1 <- make_quantmat(post_predmean_sim, 1)
#######################################################
# calculate ecdf quantile mats from each condition
# unnest post_pred_sim, using only 1 subj/dataset
unnested <-
post_predmean_sim %>%
unnest(dataset)
# calc quantiles mat for pre condition
ecdf_res_stim0 <-
unnested %>%
filter(stimulation == 0) %>%
group_by(sim_num) %>%
group_map(~ecdf(.$obs_degree )(seq(-180, 180, 1)))
stim0_ecdf_quantiles <- bind_cols(
tibble(x_val = seq(-180, 180, 1)),
as_tibble(colQuantiles(do.call(rbind, ecdf_res_stim0), probs = c(0.95, 0.5, 0.05 )))
)
# calc quantiles mat for post condition
ecdf_res_stim1 <-
unnested %>%
filter(stimulation == 1) %>%
group_by(sim_num) %>%
group_map(~ecdf(.$obs_degree )(seq(-180, 180, 1)))
stim1_ecdf_quantiles <- bind_cols(
tibble(x_val = seq(-180, 180, 1)),
as_tibble(colQuantiles(do.call(rbind, ecdf_res_stim1), probs = c(0.95, 0.5, 0.05 )))
)# blues
b_light <- "#8C9BC4"
b_light_highlight <- "#A0ADCE"
b_mid <- "#546BA9"
b_mid_highlight <- "#7385B8"
b_dark <- "#002381"
b_dark_highlight <- "#2E4B97"
#betancourt reds
r_light <- "#DCBCBC"
r_light_highlight <- "#C79999"
r_mid <- "#B97C7C"
r_mid_highlight <- "#A25050"
r_dark <- "#8F2727"
r_dark_highlight <- "#7C0000"
#######################################################
# plot histogram(density) per condition
ggplot(quantmat_cond0, aes(x = xvals)) +
geom_ribbon(aes(ymax = p0.9, ymin = p0.1), fill = r_light, alpha = 0.4) +
geom_line(aes(y = p0.5), color = r_dark, size = 1) +
geom_ribbon(data = quantmat_cond1, aes(ymax = p0.9, ymin = p0.1), fill = b_light, alpha = 0.4) +
geom_line(data = quantmat_cond1, aes(y = p0.5), color = b_dark, size = 1) +
scale_x_continuous(breaks=pretty_breaks(10)) +
labs(x = "error (degrees) [red = pre, blue = post]",
y = "count +/- quantile",
subtitle = glue("per-condition posterior pred dist (median line, 90% interval over {nrow(post_predmean_sim)} sim datasets)\n({nobs_per_cond_sim} samples/cond, per subj-level draw, per group-level mean draw)")) +
theme_cowplot()
#######################################################
# plot ecdf per condition
ggplot() +
geom_ribbon(data = stim0_ecdf_quantiles, aes(x = x_val, ymax = `95%`, ymin = `5%`), fill = "red", alpha = 0.3) +
geom_line(data = stim0_ecdf_quantiles, aes(x = x_val, y = `50%`), color = "red", size = 1) +
geom_ribbon(data = stim1_ecdf_quantiles, aes(x = x_val, ymax = `95%`, ymin = `5%`), fill = "blue", alpha = 0.3) +
geom_line(data = stim1_ecdf_quantiles, aes(x = x_val, y = `50%`), color = "blue", size = 1) +
scale_x_continuous(breaks=pretty_breaks(10)) +
geom_hline(yintercept = seq(0, 1, 0.25), linetype = "dashed", alpha = 0.2) +
labs(x = "error (degrees) [red = pre, blue = post]",
y = "cumulative prob.",
subtitle = glue("per-condition posterior pred cdf (median line, 90% interval over {nrow(post_predmean_sim)} sim datasets) \n({nobs_per_cond_sim} samples/cond, per subj-level draw, per group-level mean draw")) +
theme_cowplot()
* Notes
Predicting data, including group-level SD, shows that there is strong overlap in the pre and post predicitve dists.
Predicting data, ignoring group-level SD and predicting only the average subject, shows more of the trend of an effect.